B402A Assignment 7
Lillian Chen, Milan Filipovic, Vivian Truong
5/24/2021
Allyson Terry’s study on immunologic response to liver allograft injury
Statement of statistical hypothesis(es):
- What cytokines are significantly different between IRI- and IRI+ patients in the larger cohort of 100 liver transplant patients?
- Do these cytokines influence IRI individually?
- Do these cytokines influence IRI together?
- Are any of these cytokines significantly different at more than one timepoint?
To determine if cytokines influence IRI individually, binary logistic regression will be conducted with IRI injury as the outcome variable, with difference between timepoints in cytokine concentrations as covariates.
To address what cytokines influence IRI together, a correlation plot was constructed to visualize related cytokines.
To address whether cytokines significantly differ at more than one timepoint, separate one-way anovas with visit type as the factor variable will be conducted with each cytokine as the outcome variable. Furthermore, two-way anovas will be conducted for each cytokine concentration with visit type and IRI status as factor variables to explore if time and IRI status explain variation in cytokine concentration.
- What cytokines in the pre-transplant blood sample specifically correlate with:
- Immune cell infiltration after transplant? (Inflammation score)
- Dead cells in the liver after transplant? (Necrosis score)
Background information:
Ischemia-reperfusion injury (IRI) is when an organ loses and then regains access to blood flow, commonly occurring with organ transplantation. The immune response elicited by IRI can be severe enough to negatively affect the function and survival of the donated organ. One of two main ways that the cells in the body communicate is through cytokines, which are secreted proteins that are present in the blood of transplant patients and can influence the degree of inflammation that occurs during and after IRI. Two characteristics that doctors use to determine if a patient is IRI+ or IRI- are the inflammation score, which is the amount of immune cells in the liver after transplant and the necorsis score, which is the amount of dead cells in the liver after transplant. These characteristics of inflammation and necrosis are scored from 0 to 4, with 0 being the least severe and 4 being the most severe.
Sample data description:
The study includes 100 liver transplant patients, 42 of which are IRI+ and 58 IRI-. The analysis includes 23 cytokines. The data set also includes the inflammation score, necrosis score, and visit type for each patient.
Variable assessment:
data <- readxl::read_xlsx("./20210521_402a.xlsx", col_names = T, skip = 1)
str(data)## tibble [300 × 29] (S3: tbl_df/tbl/data.frame)
## $ Sample # : num [1:300] 1 2 3 4 5 6 7 8 9 10 ...
## $ Patient ID : num [1:300] 4 4 4 5 5 5 6 6 6 7 ...
## $ IRI Score (–,+): num [1:300] 0 0 0 0 0 0 0 0 0 0 ...
## $ Inflammation : num [1:300] 1 1 1 1 1 1 1 1 1 1 ...
## $ Necrosis : num [1:300] 0 0 0 0 0 0 0 0 0 0 ...
## $ Visit Type : chr [1:300] "PO" "LF" "W1" "PO" ...
## $ CYTO1 : chr [1:300] "3.62" "8.57" "70.98" "13.75" ...
## $ CYTO2 : chr [1:300] "18.399999999999999" "361.15" "162.97" "120.96" ...
## $ CYTO3 : num [1:300] 293 204 206 755 737 ...
## $ CYTO4 : chr [1:300] "70.069999999999993" "68.459999999999994" "12.56" "2.63" ...
## $ CYTO5 : chr [1:300] "<0.54↓" "2.76" "3.6" "7.51" ...
## $ CYTO6 : num [1:300] 101 260 1558 364 1079 ...
## $ CYTO7 : num [1:300] 22.1 3065 50.1 44.3 240.4 ...
## $ CYTO8 : chr [1:300] "<2.26↓" "3.21" "<2.26↓" "8.68" ...
## $ CYTO9 : chr [1:300] "7.41" "5.75" "20.81" "64.84" ...
## $ CYTO10 : chr [1:300] "1.99" "10.62" "9.67" "8.9700000000000006" ...
## $ CYTO11 : chr [1:300] "58.41" "187.22" "1478" "149.26" ...
## $ CYTO12 : chr [1:300] "1.56" "3.06" "1.19" "8.74" ...
## $ CYTO13 : chr [1:300] "25.57" "1210" "39.07" "127.7" ...
## $ CYTO14 : chr [1:300] "41.76" "29.2" "97.12" "169.95" ...
## $ CYTO15 : chr [1:300] "<1.45↓" "1.9" "1.98" "5.51" ...
## $ CYTO16 : chr [1:300] "<1.09↓" "1.39" "<1.09↓" "<1.09↓" ...
## $ CYTO17 : chr [1:300] "<1.38↓" "<1.38↓" "<1.38↓" "<1.38↓" ...
## $ CYTO18 : chr [1:300] "261.95" "107.07" "671.93" "2177" ...
## $ CYTO19 : chr [1:300] "5.0599999999999996" "3.49" "15.65" "32.85" ...
## $ CYTO20 : chr [1:300] "71.97" "298.99" "98.69" "203.47" ...
## $ CYTO21 : chr [1:300] "<2.07↓" "13.49" "6.1" "14.9" ...
## $ CYTO22 : num [1:300] 37.8 38.1 43.2 103.3 61.1 ...
## $ CYTO23 : num [1:300] 17.66 18.43 15.37 22.27 2.67 ...dataclean <- data %>%
rename(sample = "Sample #",
patientid = "Patient ID",
IRI = "IRI Score (–,+)",
visittype = "Visit Type") %>%
mutate(across(starts_with("CYTO"),
~ dplyr::case_when(
str_starts(.x, "<") == T ~ as.numeric(str_sub(.x, 2, -2)),
str_starts(.x, ">") == T ~ as.numeric(str_sub(.x, 2, -2)) + 1,
str_starts(.x, "1|2|3|4|5|6|7|8|9|0") == T ~ as.numeric(.x)
)))## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercion## Warning: Problem with `mutate()` input `..1`.
## ℹ NAs introduced by coercion
## ℹ Input `..1` is `across(...)`.## Warning in eval_tidy(pair$rhs, env = default_env): NAs introduced by coercionlibrary(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transposeisidata <- fread(file = "dataclean.csv", header = T, stringsAsFactors = F)
(isi_tble <- tibble(isidata))Create patient characteristics table.
#namess
library(gtsummary)
CYTO <- names(isi_tble[7:29]) #vector of cytokine variable names
demographic_table <- isi_tble %>%
mutate(IRI = ifelse( IRI == 0, "IRI-", "IRI+")) %>%
select(-patientid, -sample) %>%
tbl_summary(by = IRI,
statistic = CYTO ~
"{median} ({p25}, {p75})",
digits = all_continuous() ~ 2,
label = list(visittype ~ "Visit Type"),
missing_text = "Missing") %>%
add_overall() %>%
bold_labels() %>%
add_p() %>%
bold_p() %>%
modify_header(label ~ "**Table 1. Patient Characteristics**" )## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(CYTO)` instead of `CYTO` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.demographic_table| Table 1. Patient Characteristics | Overall, N = 3001 | IRI-, N = 1741 | IRI+, N = 1261 | p-value2 |
|---|---|---|---|---|
| Inflammation | <0.001 | |||
| 0 | 27 (9.0%) | 24 (14%) | 3 (2.4%) | |
| 1 | 177 (59%) | 126 (72%) | 51 (40%) | |
| 2 | 78 (26%) | 21 (12%) | 57 (45%) | |
| 3 | 18 (6.0%) | 3 (1.7%) | 15 (12%) | |
| Necrosis | <0.001 | |||
| 0 | 75 (25%) | 72 (41%) | 3 (2.4%) | |
| 1 | 111 (37%) | 99 (57%) | 12 (9.5%) | |
| 2 | 72 (24%) | 3 (1.7%) | 69 (55%) | |
| 3 | 36 (12%) | 0 (0%) | 36 (29%) | |
| 4 | 6 (2.0%) | 0 (0%) | 6 (4.8%) | |
| Visit Type | >0.9 | |||
| LF | 100 (33%) | 58 (33%) | 42 (33%) | |
| PO | 100 (33%) | 58 (33%) | 42 (33%) | |
| W1 | 100 (33%) | 58 (33%) | 42 (33%) | |
| CYTO1 | 8.70 (4.30, 15.69) | 8.94 (4.04, 15.62) | 8.62 (4.83, 15.55) | 0.6 |
| CYTO2 | 128.04 (66.54, 359.61) | 122.28 (66.32, 360.64) | 129.10 (69.68, 352.29) | 0.9 |
| CYTO3 | 229.10 (151.17, 324.17) | 232.78 (151.52, 326.69) | 221.18 (149.24, 312.09) | 0.7 |
| CYTO4 | 24.62 (11.54, 45.87) | 23.44 (13.08, 43.24) | 25.07 (9.51, 47.88) | 0.9 |
| CYTO5 | 3.79 (1.53, 8.85) | 3.79 (1.76, 8.32) | 3.80 (1.33, 9.34) | 0.8 |
| CYTO6 | 271.99 (152.17, 459.83) | 277.15 (172.53, 509.27) | 262.06 (140.55, 397.65) | 0.2 |
| CYTO7 | 88.35 (44.57, 243.92) | 94.28 (44.34, 276.30) | 84.80 (46.80, 198.98) | 0.7 |
| CYTO8 | 3.33 (2.01, 7.73) | 3.36 (2.00, 7.79) | 3.29 (2.10, 7.70) | 0.8 |
| CYTO9 | 15.16 (5.07, 46.98) | 15.32 (4.51, 48.28) | 14.98 (6.52, 42.99) | 0.8 |
| CYTO10 | 10.12 (6.03, 15.21) | 9.62 (5.75, 14.88) | 10.27 (6.60, 16.62) | 0.3 |
| CYTO11 | 92.06 (44.36, 178.51) | 101.12 (46.94, 189.43) | 68.12 (37.17, 158.68) | 0.037 |
| CYTO12 | 2.86 (1.58, 6.04) | 2.66 (1.58, 6.16) | 2.99 (1.56, 5.90) | 0.5 |
| CYTO13 | 149.50 (69.78, 1,568.25) | 145.42 (59.26, 1,468.50) | 168.98 (82.33, 2,155.00) | 0.3 |
| CYTO14 | 66.06 (20.62, 189.05) | 63.72 (17.14, 197.00) | 71.24 (26.99, 179.49) | 0.5 |
| CYTO15 | 3.16 (1.51, 7.11) | 3.00 (1.48, 7.26) | 3.27 (1.67, 6.22) | 0.8 |
| CYTO16 | 1.10 (0.94, 1.60) | 1.10 (0.81, 1.55) | 1.10 (0.94, 1.92) | 0.064 |
| CYTO17 | 1.38 (1.33, 1.41) | 1.39 (1.33, 1.41) | 1.34 (1.06, 1.41) | 0.018 |
| CYTO18 | 431.17 (131.20, 1,351.75) | 388.58 (109.13, 1,363.25) | 491.25 (189.98, 1,316.75) | 0.3 |
| CYTO19 | 8.79 (3.79, 24.11) | 7.98 (3.34, 25.27) | 10.02 (4.52, 21.50) | 0.5 |
| CYTO20 | 152.46 (68.05, 306.70) | 137.80 (61.75, 298.48) | 163.06 (70.80, 336.67) | 0.3 |
| CYTO21 | 5.24 (2.79, 12.36) | 5.24 (2.59, 13.66) | 5.26 (3.28, 11.09) | 0.8 |
| CYTO22 | 58.27 (35.16, 107.86) | 49.77 (31.09, 92.57) | 75.42 (37.61, 117.85) | 0.005 |
| CYTO23 | 20.20 (13.20, 30.07) | 18.40 (12.97, 28.57) | 21.86 (13.45, 32.78) | 0.065 |
|
1
n (%); Median (IQR)
2
Pearson's Chi-squared test; Fisher's exact test; Wilcoxon rank sum test
|
||||
It should be noted that a significant difference between IRI- and IRI+ patients was detected for necrosis score (p <.001).
The significance test for cytokine concentration was conducted using a nonparametric method, so it is more robust to assumptions of normality. CYTO11, CYTO17, and CYTO22 were found to be significantly different. Next we visualize the distribution of cytokine concentration.
isi_tble %>%
select(CYTO1:CYTO6) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO7:CYTO12) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO13:CYTO18) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO19:CYTO23) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## my attempt to facet in one plot. Does not look good
isi_tble %>% select(starts_with("CYT")) %>%
gather() %>%
ggplot(aes(value)) +
facet_wrap(~ key, scales ="free") +
# facet_grid_paginate(~ key, ncol = 3, nrow = 3, page = 2)
geom_boxplot()
The Cytokine distributions are characterized by heavy right skew.
Now log transform all Cytokines
isi_tble %>%
select(CYTO1:CYTO6) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram() +
labs(title = "Log Transformed") +
xlab("log pg/mL")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO7:CYTO12) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram() +
labs(title = "Log Transformed") +
xlab("log pg/mL")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO13:CYTO18) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram() +
labs(title = "Log Transformed") +
xlab("log pg/mL")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO19:CYTO23) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_histogram() +
labs(title = "Log Transformed") +
xlab("log pg/mL")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
isi_tble %>%
select(CYTO1:CYTO6) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_boxplot() +
labs(title = "Log Transformed") +
xlab("log pg/mL")
isi_tble %>%
select(CYTO7:CYTO12) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_boxplot() +
labs(title = "Log Transformed") +
xlab("log pg/mL")
isi_tble %>%
select(CYTO13:CYTO18) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_boxplot() +
labs(title = "Log Transformed") +
xlab("log pg/mL")
isi_tble %>%
select(CYTO19:CYTO23) %>%
gather() %>%
ggplot(aes(log(value))) +
facet_wrap(~ key, nrow = 3, scales = "free") +
geom_boxplot() +
labs(title = "Log Transformed") +
xlab("log pg/mL")
- Among the log transformed cytokine values, it is seen that that overall normality is improved.
- CYTO12 and CYTO 19 appears to have a similar distribution when looking at the histogram, with large left skew and lack of density left of the peak density.
- CYTO13 appears as if it may have another peak in the right tail. It is right skewed.
- CYTO 15 appears to have several outliers in the higher range.
- CYTO16 appears right skewed.
- CYTO18 appears to have lower outliers.
- CYTO17 has many values below level of detection, therefore the transformation did not affect overall normality.
- CYTO21 appears to have a large right skew, though it should be noted that many values were below level of detection in CYTO 21.
It should be noted that due to many observations below detection, CytO17 should not be included in the linear models that follow.
Cytokines Trending With Time by Anova
First look at how cytokine values trend overtime
library(ggpubr)
#Reordering the visittype
isi_tble$visittype <- factor(isi_tble$visittype, levels = c("PO", "LF", "W1"))
#One-way anova at each timepoint
for(i in c(7:29)){
p <- isi_tble %>%
ggplot() +
geom_boxplot(aes(x = visittype, y = log(isi_tble[[i]]), color = visittype)) +
stat_compare_means(aes( x = visittype, y = log(isi_tble[[i]])), method = "anova",
label.y = (max(log(isi_tble[[i]])) + 1)) + #positions the anova stat
# stat_compare_means(aes( x = visittype, y = log(isi_tble[[i]])), label = "p.signif", method = "t.test",
# ref.group = ".all.") +
geom_hline(yintercept = mean(log(isi_tble[[i]])), linetype = 2) + #adds horizontal line for mean
ggtitle(paste("Histogram of Log", names(isi_tble)[i], "by Visit Type")) +
ylab("pg/mL") +
xlab("Visit Type") +
theme_classic()
print(p)
}## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
## Warning: Use of `isi_tble[[i]]` is discouraged. Use `.data[[i]]` instead.
Only CYTO4, CYTO16, and CYTO17 concentrations are not significantly different at the different visits by one way anova. In all other cytokines, an overall difference in mean concentration was dtected. The overall mean at all timepoints is displayed with a horizontal line.
CYTO2 at the LF timepoint appears significantly larger. CYTO7 at the LF timepoint appears significantly larger. CYTO13 at the LF timepoint appears significantly larger. In many of the cytokines, the concentration at the LF timepoint is larger, suggesting an elevated immune response during the LF timepoint. These cytokines include CYTO2, CYTO6, CYTO 7, CYTO13, CYTO20. The concentration of CYTO22 appears to decrease at each timepoint.
To investigate this further, pairwise comparisons will be conducted using Tukey’s HSD.
# Log Transforming Data to use in anova modelling.
liri_tble <- isi_tble %>% select(starts_with("CYT")) %>%
log() %>% #log transforms the cytokine data frame
mutate(visittype = isi_tble$visittype, IRI = isi_tble$IRI, ID = isi_tble$patientid) #add the factor variables back in
responseList <- names(liri_tble)[c(1:23)]
modelList <- lapply(responseList, function(resp) {
mF <- formula(paste(resp, " ~ visittype"))
aov(mF, data = liri_tble)
})
HSDList <- lapply(modelList, function(mod) {
TukeyHSD(mod)
})
HSDList## [[1]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.06840745 -0.2558110 0.39262592 0.8728010
## W1-PO -0.38474438 -0.7089628 -0.06052591 0.0152141
## W1-LF -0.45315183 -0.7773703 -0.12893337 0.0031873
##
##
## [[2]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 1.4051474 1.1024736 1.7078211 0.00e+00
## W1-PO -0.5458561 -0.8485298 -0.2431823 8.55e-05
## W1-LF -1.9510034 -2.2536772 -1.6483296 0.00e+00
##
##
## [[3]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.34268195 -0.5428843 -0.1424796 0.0002074
## W1-PO -0.38544874 -0.5856511 -0.1852464 0.0000249
## W1-LF -0.04276679 -0.2429691 0.1574355 0.8698349
##
##
## [[4]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.333376411 -0.08740419 0.75415701 0.1503745
## W1-PO -0.007703802 -0.42848440 0.41307680 0.9989752
## W1-LF -0.341080213 -0.76186081 0.07970039 0.1377927
##
##
## [[5]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.3420904 -0.9707344 0.2865535 0.4065536
## W1-PO -0.8624686 -1.4911126 -0.2338247 0.0039032
## W1-LF -0.5203782 -1.1490221 0.1082658 0.1267213
##
##
## [[6]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.79897356 0.5359254 1.0620217 0.0000000
## W1-PO 0.06966596 -0.1933822 0.3327141 0.8071605
## W1-LF -0.72930760 -0.9923558 -0.4662594 0.0000000
##
##
## [[7]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 1.6747311 1.31697119 2.032491 0.0000000
## W1-PO 0.2926001 -0.06515987 0.650360 0.1330189
## W1-LF -1.3821311 -1.73989102 -1.024371 0.0000000
##
##
## [[8]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.1599651 -0.5705408 0.2506106 0.6295108
## W1-PO -0.9467543 -1.3573300 -0.5361786 0.0000003
## W1-LF -0.7867892 -1.1973649 -0.3762136 0.0000273
##
##
## [[9]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.1335385 -0.3132586 0.5803357 0.7613020
## W1-PO -0.8037803 -1.2505775 -0.3569831 0.0000894
## W1-LF -0.9373188 -1.3841160 -0.4905217 0.0000039
##
##
## [[10]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.5342223 0.2846603 0.7837843 0.0000024
## W1-PO 0.1231167 -0.1264454 0.3726787 0.4768492
## W1-LF -0.4111056 -0.6606677 -0.1615436 0.0003772
##
##
## [[11]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.6541453 -1.0230812 -0.2852093 0.0001152
## W1-PO -0.5356416 -0.9045775 -0.1667056 0.0020579
## W1-LF 0.1185037 -0.2504323 0.4874397 0.7298669
##
##
## [[12]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.4644485 -0.9302121 0.001315079 0.0508349
## W1-PO -0.8359971 -1.3017606 -0.370233482 0.0000930
## W1-LF -0.3715486 -0.8373121 0.094215016 0.1465507
##
##
## [[13]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 3.6265834 3.2519318 4.00123497 0.0000000
## W1-PO -0.3992638 -0.7739155 -0.02461223 0.0336001
## W1-LF -4.0258472 -4.4004988 -3.65119559 0.0000000
##
##
## [[14]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.4232374 -0.2033776 1.04985235 0.2510970
## W1-PO -0.3004576 -0.9270725 0.32615740 0.4966867
## W1-LF -0.7236950 -1.3503099 -0.09707998 0.0188714
##
##
## [[15]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.2633009 -0.02313782 0.5497397 0.0789346
## W1-PO -0.6853100 -0.97174879 -0.3988713 0.0000001
## W1-LF -0.9486110 -1.23504972 -0.6621722 0.0000000
##
##
## [[16]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.17846371 -0.424894948 0.06796753 0.2046683
## W1-PO 0.07083474 -0.175596497 0.31726598 0.7770205
## W1-LF 0.24929845 0.002867211 0.49572969 0.0466803
##
##
## [[17]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.005186610 -0.10270436 0.09233114 0.9913845
## W1-PO -0.004979775 -0.10249752 0.09253797 0.9920552
## W1-LF 0.000206835 -0.09731091 0.09772458 0.9999862
##
##
## [[18]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.0343640 -0.698676 0.76740397 0.9933003
## W1-PO -0.7552324 -1.488272 -0.02219241 0.0417514
## W1-LF -0.7895964 -1.522636 -0.05655641 0.0312749
##
##
## [[19]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.09632616 -0.3560285 0.5486808 0.8705987
## W1-PO -0.53718522 -0.9895398 -0.0848306 0.0151257
## W1-LF -0.63351138 -1.0858660 -0.1811568 0.0031171
##
##
## [[20]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.6492907 0.2646393 1.0339421 0.0002590
## W1-PO -0.6123829 -0.9970343 -0.2277315 0.0006208
## W1-LF -1.2616736 -1.6463250 -0.8770222 0.0000000
##
##
## [[21]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO 0.04883302 -0.2383852 0.3360513 0.9154208
## W1-PO -0.77089036 -1.0581086 -0.4836721 0.0000000
## W1-LF -0.81972338 -1.1069416 -0.5325051 0.0000000
##
##
## [[22]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.1014427 -0.3726096 0.1697241 0.6525960
## W1-PO -0.6548417 -0.9260086 -0.3836748 0.0000001
## W1-LF -0.5533990 -0.8245658 -0.2822321 0.0000073
##
##
## [[23]]
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = mF, data = liri_tble)
##
## $visittype
## diff lwr upr p adj
## LF-PO -0.6466722 -0.82377924 -0.4695653 0.0000000
## W1-PO 0.1462652 -0.03084183 0.3233722 0.1279238
## W1-LF 0.7929374 0.61583041 0.9700444 0.0000000The Cytokines in which concentration at timepoint PO is significantly different from LF and W1, but LF and W1 are NOT significantly different from each other: CYTO6, CYTO 7, CYTO10, CYTO13, and CYTO23. It should be noted that in CYTO23, concentration at LF was actually significantly lower than LF and W1. With the exception of CYTO23, these cytokines suggest that the LF timepoint is associated with an increased immune response.
CYTO 2, CYTO13 and CYTO20 were significantly different for each pairwise comparison.
To explore how cytokine concentrations trend over time in IRI positive and negative patients, interaction plots between IRI status and visit are visualized below. Here median values are plotted over the visit types. It should be noted that CYTO4, CYTO16, and CYTO17 were not included because significance was not found in the anova test.
#need to add titles here
#I use median instead of mean due to outliers in data set
for(i in c(1:20)){
iri_tble2<- liri_tble %>% select(starts_with("CYTO"), -CYTO4, -CYTO16, -CYTO17)
interaction.plot(x.factor = liri_tble$visittype,
response = liri_tble[[i]],
trace.factor = liri_tble$IRI,
fun = median,
ylab = "log(pg/mL)",
xlab = paste(names(iri_tble2)[i], "Visit Type"),
col = c("red", "blue"),
trace.label = "IRI Status")
}
The interaction plots indicate possible interaction effect between IRI status and visit type in CYTO5, 7, 8, 14, 19 and 21.
To explore whether IRI status explains variation in cytokine concentrations, two-way anovas will be conducted with IRI status and visit type as factor levels. A separate anova for each cytokine concentration will be conducted.
Running 2-Way ANOVA
twoway_modelList <- lapply(responseList, function(resp) {
mF <- formula(paste(resp, " ~ visittype*IRI"))
anova(lm(mF, data = isi_tble))
})
twoway_modelList## [[1]]
## Analysis of Variance Table
##
## Response: CYTO1
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 10755178 5377589 0.9785 0.3771
## IRI 1 7664312 7664312 1.3946 0.2386
## visittype:IRI 2 15242052 7621026 1.3867 0.2515
## Residuals 294 1615708005 5495605
##
## [[2]]
## Analysis of Variance Table
##
## Response: CYTO2
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 15252800 7626400 51.8650 <2e-16 ***
## IRI 1 126116 126116 0.8577 0.3551
## visittype:IRI 2 149221 74611 0.5074 0.6026
## Residuals 294 43230701 147043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[3]]
## Analysis of Variance Table
##
## Response: CYTO3
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 855087 427543 13.7892 1.888e-06 ***
## IRI 1 1 1 0.0000 0.9946
## visittype:IRI 2 1846 923 0.0298 0.9707
## Residuals 294 9115680 31006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[4]]
## Analysis of Variance Table
##
## Response: CYTO4
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 1603 801.51 0.3167 0.7288
## IRI 1 346 345.51 0.1365 0.7120
## visittype:IRI 2 809 404.36 0.1598 0.8524
## Residuals 294 744146 2531.11
##
## [[5]]
## Analysis of Variance Table
##
## Response: CYTO5
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 25249 12624.5 1.9916 0.1383
## IRI 1 11428 11427.6 1.8028 0.1804
## visittype:IRI 2 11342 5670.9 0.8946 0.4099
## Residuals 294 1863628 6338.9
##
## [[6]]
## Analysis of Variance Table
##
## Response: CYTO6
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 3706755 1853377 12.1491 8.52e-06 ***
## IRI 1 43557 43557 0.2855 0.5935
## visittype:IRI 2 63807 31904 0.2091 0.8114
## Residuals 294 44850392 152552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[7]]
## Analysis of Variance Table
##
## Response: CYTO7
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 12968524 6484262 30.1170 1.263e-12 ***
## IRI 1 27613 27613 0.1283 0.7205
## visittype:IRI 2 180210 90105 0.4185 0.6584
## Residuals 294 63298826 215302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[8]]
## Analysis of Variance Table
##
## Response: CYTO8
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 27674 13837.0 2.3348 0.09862 .
## IRI 1 2256 2255.6 0.3806 0.53776
## visittype:IRI 2 6838 3418.8 0.5769 0.56228
## Residuals 294 1742357 5926.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[9]]
## Analysis of Variance Table
##
## Response: CYTO9
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 50490 25245.1 9.2506 0.000127 ***
## IRI 1 691 691.0 0.2532 0.615195
## visittype:IRI 2 54 27.1 0.0099 0.990113
## Residuals 294 802337 2729.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[10]]
## Analysis of Variance Table
##
## Response: CYTO10
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 4059 2029.30 5.2183 0.00593 **
## IRI 1 0 0.09 0.0002 0.98755
## visittype:IRI 2 140 69.80 0.1795 0.83579
## Residuals 294 114332 388.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[11]]
## Analysis of Variance Table
##
## Response: CYTO11
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 491474 245737 3.2097 0.04179 *
## IRI 1 193697 193697 2.5300 0.11278
## visittype:IRI 2 24977 12489 0.1631 0.84957
## Residuals 294 22509036 76561
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[12]]
## Analysis of Variance Table
##
## Response: CYTO12
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 4741 2370.51 9.4276 0.0001075 ***
## IRI 1 268 268.18 1.0666 0.3025683
## visittype:IRI 2 54 27.04 0.1075 0.8980720
## Residuals 294 73924 251.44
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[13]]
## Analysis of Variance Table
##
## Response: CYTO13
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 4987383986 2493691993 87.1907 <2e-16 ***
## IRI 1 399713 399713 0.0140 0.906
## visittype:IRI 2 3235162 1617581 0.0566 0.945
## Residuals 294 8408523863 28600421
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[14]]
## Analysis of Variance Table
##
## Response: CYTO14
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 112269 56134 0.0969 0.9077
## IRI 1 761089 761089 1.3141 0.2526
## visittype:IRI 2 1982831 991415 1.7118 0.1823
## Residuals 294 170270008 579150
##
## [[15]]
## Analysis of Variance Table
##
## Response: CYTO15
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 1757 878.58 4.8096 0.008804 **
## IRI 1 103 103.03 0.5640 0.453256
## visittype:IRI 2 181 90.54 0.4956 0.609694
## Residuals 294 53706 182.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[16]]
## Analysis of Variance Table
##
## Response: CYTO16
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 108.8 54.399 1.5242 0.2195
## IRI 1 49.6 49.570 1.3889 0.2395
## visittype:IRI 2 53.0 26.523 0.7432 0.4765
## Residuals 294 10492.5 35.689
##
## [[17]]
## Analysis of Variance Table
##
## Response: CYTO17
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 0.028 0.01378 0.0605 0.94129
## IRI 1 1.343 1.34292 5.8991 0.01575 *
## visittype:IRI 2 0.020 0.00998 0.0438 0.95713
## Residuals 294 66.928 0.22765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[18]]
## Analysis of Variance Table
##
## Response: CYTO18
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 53859060 26929530 3.3775 0.03546 *
## IRI 1 235729 235729 0.0296 0.86360
## visittype:IRI 2 21422916 10711458 1.3434 0.26254
## Residuals 294 2344107142 7973154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[19]]
## Analysis of Variance Table
##
## Response: CYTO19
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 6854 3427.2 5.6209 0.004022 **
## IRI 1 31 31.4 0.0515 0.820546
## visittype:IRI 2 336 167.8 0.2752 0.759641
## Residuals 294 179262 609.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[20]]
## Analysis of Variance Table
##
## Response: CYTO20
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 4450502 2225251 7.6730 0.0005646 ***
## IRI 1 242530 242530 0.8363 0.3612118
## visittype:IRI 2 550547 275273 0.9492 0.3882387
## Residuals 294 85262873 290010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[21]]
## Analysis of Variance Table
##
## Response: CYTO21
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 5746 2872.87 13.1263 3.464e-06 ***
## IRI 1 85 84.98 0.3883 0.5337
## visittype:IRI 2 18 9.19 0.0420 0.9589
## Residuals 294 64346 218.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[22]]
## Analysis of Variance Table
##
## Response: CYTO22
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 264059 132030 15.7549 3.163e-07 ***
## IRI 1 32521 32521 3.8807 0.04978 *
## visittype:IRI 2 7358 3679 0.4390 0.64510
## Residuals 294 2463786 8380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## [[23]]
## Analysis of Variance Table
##
## Response: CYTO23
## Df Sum Sq Mean Sq F value Pr(>F)
## visittype 2 13467 6733.6 45.7942 < 2e-16 ***
## IRI 1 584 583.6 3.9687 0.04728 *
## visittype:IRI 2 263 131.7 0.8954 0.40955
## Residuals 294 43230 147.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A significant interaction is not observed in any of the models between IRI and visittype. IRI effect is significant explaining the variation in CYTO23 and CYTO22, which indicates IRI status is associated with CYTO23 and CYTO22 given visit type.
Running Logistic Regression Models - Main Effects and Select Interactions
Correlation Plots
library(ggcorrplot)
corr.main <- round(cor(logdataclean[,c(3, 7:29)]), 3)
corr.change <- round(cor(changesubset[,3:26]), 3)
pcorr.change <- round(cor_pmat(changesubset[,3:26]), 3)
corr.diff1 <- round(cor(diff1subset[,3:26]), 3)
pcorr.diff1 <- round(cor_pmat(diff1subset[,3:26]), 3)
corr.diff2 <- round(cor(diff2subset[,3:26]), 3)
pcorr.diff2 <- round(cor_pmat(diff2subset[,3:26]), 3)
ggcorrplot(corr.main, outline.col = "white", type = "lower")
ggcorrplot(corr.change, outline.col = "white", type = "lower")
ggcorrplot(corr.diff1, outline.col = "white", type = "lower")
ggcorrplot(corr.diff2, outline.col = "white", type = "lower")
# did not include interaction model for change models due to not all models converging
main.table | Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.5 | -5.2, 2.1 | 0.4 |
| CYTO1 | 0.41 | 0.04, 0.85 | 0.054 |
| CYTO16 | 0.30 | -0.03, 0.65 | 0.081 |
| CYTO18 | 0.23 | 0.03, 0.50 | 0.046 |
| CYTO22 | 0.41 | 0.07, 0.77 | 0.020 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
po.table | Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.4 | -8.6, 5.6 | 0.7 |
| CYTO3 | -0.54 | -1.2, 0.07 | 0.091 |
| CYTO15 | -0.92 | -2.0, 0.03 | 0.070 |
| CYTO16 | 0.67 | -0.08, 1.5 | 0.094 |
| CYTO23 | 0.76 | -0.05, 1.7 | 0.077 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
lf.table | Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.9 | -15, 11 | 0.8 |
| CYTO1 | 1.2 | 0.08, 2.5 | 0.044 |
| CYTO2 | -1.2 | -2.2, -0.28 | 0.017 |
| CYTO11 | -0.55 | -1.1, -0.03 | 0.049 |
| CYTO13 | 0.45 | -0.02, 1.0 | 0.075 |
| CYTO14 | -1.2 | -2.7, -0.18 | 0.054 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
w1.table | Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -0.18 | -11, 10 | >0.9 |
| CYTO9 | -1.1 | -2.2, -0.06 | 0.048 |
| CYTO13 | -0.58 | -1.3, 0.06 | 0.088 |
| CYTO18 | 0.73 | 0.03, 1.6 | 0.070 |
| CYTO22 | 1.6 | 0.35, 2.9 | 0.016 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
stepint.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.7 | -2.9, -0.55 | 0.004 |
| CYTO16 | 0.22 | -0.01, 0.45 | 0.063 |
| CYTO22 | 0.22 | 0.02, 0.42 | 0.031 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
steppoint.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -7.6 | -16, -0.62 | 0.045 |
| CYTO3 | 0.80 | -0.15, 1.9 | 0.11 |
| CYTO15 | 1.9 | -0.66, 5.0 | 0.2 |
| CYTO23 | 0.23 | -0.70, 1.2 | 0.6 |
| CYTO3 * CYTO15 | -0.45 | -0.90, -0.05 | 0.036 |
| CYTO15 * CYTO23 | 0.35 | -0.08, 0.82 | 0.13 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
steplfint.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 10 | 0.50, 22 | 0.051 |
| CYTO1 | 0.59 | -0.14, 1.4 | 0.12 |
| CYTO2 | -0.83 | -1.5, -0.27 | 0.006 |
| CYTO11 | -1.5 | -3.2, -0.08 | 0.050 |
| CYTO13 | -0.12 | -1.0, 0.63 | 0.8 |
| CYTO14 | -0.27 | -0.62, 0.07 | 0.12 |
| CYTO11 * CYTO13 | 0.09 | -0.02, 0.23 | 0.13 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
stepw1int.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -17 | -36, -2.1 | 0.045 |
| CYTO9 | -1.3 | -2.2, -0.57 | 0.001 |
| CYTO13 | 1.4 | -1.2, 4.4 | 0.3 |
| CYTO18 | -1.6 | -3.4, 0.06 | 0.066 |
| CYTO22 | 6.4 | 2.2, 12 | 0.008 |
| CYTO13 * CYTO18 | 0.36 | 0.10, 0.67 | 0.011 |
| CYTO13 * CYTO22 | -0.81 | -1.6, -0.15 | 0.030 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
changemain.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -0.25 | -0.67, 0.17 | 0.2 |
| CYTO9.diffs | -0.44 | -1.0, 0.04 | 0.084 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
changemainbase.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.5 | -7.4, 4.3 | 0.6 |
| CYTO9.diffs | -0.55 | -1.2, 0.08 | 0.10 |
| CYTO1 | 0.81 | -0.01, 1.7 | 0.056 |
| CYTO3 | -0.54 | -1.1, -0.05 | 0.035 |
| CYTO15 | -0.82 | -1.6, -0.06 | 0.039 |
| CYTO16 | 0.71 | 0.09, 1.4 | 0.031 |
| CYTO23 | 0.78 | 0.14, 1.5 | 0.020 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
changemain1.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.8 | -3.6, -0.10 | 0.045 |
| CYTO2.diffs | -0.57 | -1.2, -0.03 | 0.051 |
| CYTO3.diffs | 0.64 | -0.04, 1.4 | 0.076 |
| CYTO4.diffs | 0.51 | 0.06, 1.0 | 0.035 |
| CYTO9.diffs | 1.4 | 0.29, 2.8 | 0.024 |
| CYTO10.diffs | 0.62 | 0.00, 1.4 | 0.068 |
| CYTO14.diffs | -0.86 | -1.7, -0.18 | 0.022 |
| CYTO19.diffs | -1.2 | -2.5, -0.12 | 0.060 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
changemain2.table| Characteristic | log(OR)1 | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | -1.5 | -3.8, 0.74 | 0.2 |
| CYTO2.diffs | 0.44 | 0.01, 0.92 | 0.058 |
| CYTO9.diffs | -1.4 | -2.3, -0.61 | 0.001 |
|
1
OR = Odds Ratio, CI = Confidence Interval
|
|||
Trying to run logistic model with patient ID as a random effect.
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod carlibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:gtsummary':
##
## select## The following object is masked from 'package:dplyr':
##
## select#liri_tble[ ,c(1:23)] <- scale(liri_tble[ ,c(1:23)])
CYTO <- names(liri_tble[1:23])
eqn <- as.formula(paste("IRI ~ (1 | ID)", paste(CYTO, collapse = " + "), sep = "+"))
liri_tble$visittype <- as.numeric(liri_tble$visittype)
#logistic mixed model using time and intercept as a mrandom effect for each patient ID
random <- glmer(eqn, data = liri_tble, family = binomial, nAGQ = 0, control=glmerControl(optimizer="bobyqa"))
summary(random)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
## Family: binomial ( logit )
## Formula: IRI ~ (1 | ID) + CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 +
## CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 +
## CYTO14 + CYTO15 + CYTO16 + CYTO17 + CYTO18 + CYTO19 + CYTO20 +
## CYTO21 + CYTO22 + CYTO23
## Data: liri_tble
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 230.1 322.7 -90.1 180.1 275
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.20474 -0.11081 -0.07526 0.11669 0.21861
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 95.95 9.795
## Number of obs: 300, groups: ID, 100
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.32589 10.40342 -0.127 0.899
## CYTO1 0.75603 1.81737 0.416 0.677
## CYTO2 0.22492 1.21913 0.184 0.854
## CYTO3 -0.16697 1.34745 -0.124 0.901
## CYTO4 -0.09345 0.72447 -0.129 0.897
## CYTO5 -0.19017 0.55815 -0.341 0.733
## CYTO6 -0.35968 1.01224 -0.355 0.722
## CYTO7 0.06613 0.66253 0.100 0.920
## CYTO8 -0.01857 0.99660 -0.019 0.985
## CYTO9 -0.36080 1.85374 -0.195 0.846
## CYTO10 -0.02828 1.34454 -0.021 0.983
## CYTO11 -0.30705 0.92635 -0.331 0.740
## CYTO12 0.24520 0.84513 0.290 0.772
## CYTO13 0.08171 0.62510 0.131 0.896
## CYTO14 -0.41407 1.14959 -0.360 0.719
## CYTO15 -0.44441 1.62660 -0.273 0.785
## CYTO16 0.46918 1.39166 0.337 0.736
## CYTO17 -2.20861 3.94256 -0.560 0.575
## CYTO18 0.29932 0.73467 0.407 0.684
## CYTO19 0.09505 1.50139 0.063 0.950
## CYTO20 0.03489 1.19422 0.029 0.977
## CYTO21 -0.28314 1.51606 -0.187 0.852
## CYTO22 0.64339 1.31607 0.489 0.625
## CYTO23 0.34889 1.68086 0.208 0.836##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnalysis of data set:
Question 2: What cytokines in the pre-transplant blood sample specifically correlate with:
- Immune cell infiltration after transplant? (Inflammation score)
Since we are looking at cytokines in pre-transplant blood samples, we first filter for pre-operation data.
# filter for preoperation data
dataPO <- dataclean %>%
filter(visittype == "PO")
# crosstabs on preoperation data
xtabs(~ Inflammation, data = dataPO)## Inflammation
## 0 1 2 3
## 9 59 26 6xtabs(~ Necrosis, data = dataPO)## Necrosis
## 0 1 2 3 4
## 25 37 24 12 2We reassign the numerical inflammation scores to reflect negative and positive scores.
dataPO1 <- dataPO
# reassigning inflammation scores 0=negative and 2=positive
dataPO1$Inflammation[dataPO$Inflammation == 0] <- 0
dataPO1$Inflammation[dataPO$Inflammation == 1] <- 0
dataPO1$Inflammation[dataPO$Inflammation == 2] <- 2
dataPO1$Inflammation[dataPO$Inflammation == 3] <- 2
dataPO1$Inflammation[dataPO$Inflammation == 4] <- 2
# reassigning necrosis scores 0=negative and 2=positive
dataPO1$Necrosis[dataPO$Necrosis == 0] <- 0
dataPO1$Necrosis[dataPO$Necrosis == 1] <- 0
dataPO1$Necrosis[dataPO$Necrosis == 2] <- 2
dataPO1$Necrosis[dataPO$Necrosis == 3] <- 2
dataPO1$Necrosis[dataPO$Necrosis == 4] <- 2
# crosstabs to check that reassigned values match original data
xtabs(~ Inflammation, data = dataPO1)## Inflammation
## 0 2
## 68 32xtabs(~ Necrosis, data = dataPO1)## Necrosis
## 0 2
## 62 38We attempt an ordinal logistic regression model with the numeric inflammation scores. The confidence interval for cytokine 23 does not include 0, which suggests that it may be statistically significant.
library(MASS)
# releveled inflammation as factor
dataPO$Inflammation <- factor(dataPO$Inflammation, levels = c("0", "1", "2", "3", "4"))
dataPO$Inflammation## [1] 1 1 1 1 2 2 2 1 3 1 1 1 3 0 2 3 2 2 1 1 2 1 2 1 2 1 1 1 2 1 2 2 1 1 1 1 1
## [38] 1 1 1 3 1 2 1 1 1 2 1 3 1 1 2 2 1 1 1 3 1 1 1 2 1 2 1 2 2 1 2 1 1 0 1 1 1
## [75] 1 1 1 2 1 1 0 1 2 1 1 1 0 2 1 2 1 1 0 2 0 0 0 1 0 1
## Levels: 0 1 2 3 4# ordinal logistic regression model
m1 <- polr(Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess=TRUE)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(m1)## Call:
## polr(formula = Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 +
## CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 +
## CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 +
## CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## CYTO1 0.0990972 0.0607496 1.63124
## CYTO2 -0.0041909 0.0037136 -1.12854
## CYTO3 -0.0015958 0.0012598 -1.26675
## CYTO4 -0.0185309 0.0103424 -1.79174
## CYTO5 0.0010193 0.0018938 0.53824
## CYTO6 0.0011034 0.0011182 0.98678
## CYTO7 0.0012252 0.0012310 0.99532
## CYTO8 0.0007049 0.0140244 0.05026
## CYTO9 0.0081579 0.0136709 0.59674
## CYTO10 -0.0244423 0.0350566 -0.69722
## CYTO11 -0.0026796 0.0015254 -1.75668
## CYTO12 -0.0264420 0.0223970 -1.18060
## CYTO13 -0.0005690 0.0007607 -0.74808
## CYTO14 -0.0046667 0.0035552 -1.31264
## CYTO15 0.0227108 0.0506707 0.44820
## CYTO16 0.2700851 0.1510910 1.78757
## CYTO18 -0.0003892 0.0002172 -1.79196
## CYTO19 0.0405524 0.0356837 1.13644
## CYTO20 -0.0008111 0.0007264 -1.11657
## CYTO21 0.0233837 0.0277171 0.84366
## CYTO22 -0.0004304 0.0024876 -0.17303
## CYTO23 0.0578270 0.0223748 2.58447
##
## Intercepts:
## Value Std. Error t value
## 0|1 -2.0028 0.6858 -2.9204
## 1|2 1.8115 0.6294 2.8781
## 2|3 4.4719 0.8403 5.3217
## 3|4 23.4397 0.8402 27.8988
##
## Residual Deviance: 170.8675
## AIC: 222.8675## store table
(ctable1 <- coef(summary(m1)))## Value Std. Error t value
## CYTO1 0.0990972489 0.0607495853 1.63124157
## CYTO2 -0.0041909214 0.0037135846 -1.12853801
## CYTO3 -0.0015957931 0.0012597556 -1.26674817
## CYTO4 -0.0185309221 0.0103423908 -1.79174453
## CYTO5 0.0010193281 0.0018938064 0.53824307
## CYTO6 0.0011033743 0.0011181528 0.98678309
## CYTO7 0.0012252488 0.0012310094 0.99532041
## CYTO8 0.0007049127 0.0140243958 0.05026332
## CYTO9 0.0081579103 0.0136709009 0.59673538
## CYTO10 -0.0244423105 0.0350565760 -0.69722469
## CYTO11 -0.0026795576 0.0015253564 -1.75667638
## CYTO12 -0.0264420415 0.0223970444 -1.18060406
## CYTO13 -0.0005690480 0.0007606737 -0.74808429
## CYTO14 -0.0046667285 0.0035552334 -1.31263633
## CYTO15 0.0227107980 0.0506706904 0.44820384
## CYTO16 0.2700850984 0.1510910309 1.78756540
## CYTO18 -0.0003892376 0.0002172136 -1.79195821
## CYTO19 0.0405524257 0.0356837272 1.13644030
## CYTO20 -0.0008110602 0.0007263833 -1.11657342
## CYTO21 0.0233837014 0.0277170687 0.84365709
## CYTO22 -0.0004304278 0.0024876217 -0.17302785
## CYTO23 0.0578269674 0.0223747943 2.58446923
## 0|1 -2.0028022044 0.6857895393 -2.92043271
## 1|2 1.8114728183 0.6293997691 2.87809578
## 2|3 4.4719234340 0.8403158047 5.32171763
## 3|4 23.4396649994 0.8401688913 27.89875374## calculate and store p values
p1 <- pnorm(abs(ctable1[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable1 <- cbind(ctable1, "p value" = p1))## Value Std. Error t value p value
## CYTO1 0.0990972489 0.0607495853 1.63124157 1.028394e-01
## CYTO2 -0.0041909214 0.0037135846 -1.12853801 2.590928e-01
## CYTO3 -0.0015957931 0.0012597556 -1.26674817 2.052453e-01
## CYTO4 -0.0185309221 0.0103423908 -1.79174453 7.317390e-02
## CYTO5 0.0010193281 0.0018938064 0.53824307 5.904093e-01
## CYTO6 0.0011033743 0.0011181528 0.98678309 3.237490e-01
## CYTO7 0.0012252488 0.0012310094 0.99532041 3.195805e-01
## CYTO8 0.0007049127 0.0140243958 0.05026332 9.599126e-01
## CYTO9 0.0081579103 0.0136709009 0.59673538 5.506841e-01
## CYTO10 -0.0244423105 0.0350565760 -0.69722469 4.856622e-01
## CYTO11 -0.0026795576 0.0015253564 -1.75667638 7.897299e-02
## CYTO12 -0.0264420415 0.0223970444 -1.18060406 2.377601e-01
## CYTO13 -0.0005690480 0.0007606737 -0.74808429 4.544093e-01
## CYTO14 -0.0046667285 0.0035552334 -1.31263633 1.893055e-01
## CYTO15 0.0227107980 0.0506706904 0.44820384 6.540061e-01
## CYTO16 0.2700850984 0.1510910309 1.78756540 7.384615e-02
## CYTO18 -0.0003892376 0.0002172136 -1.79195821 7.313966e-02
## CYTO19 0.0405524257 0.0356837272 1.13644030 2.557723e-01
## CYTO20 -0.0008110602 0.0007263833 -1.11657342 2.641768e-01
## CYTO21 0.0233837014 0.0277170687 0.84365709 3.988611e-01
## CYTO22 -0.0004304278 0.0024876217 -0.17302785 8.626295e-01
## CYTO23 0.0578269674 0.0223747943 2.58446923 9.752901e-03
## 0|1 -2.0028022044 0.6857895393 -2.92043271 3.495457e-03
## 1|2 1.8114728183 0.6293997691 2.87809578 4.000836e-03
## 2|3 4.4719234340 0.8403158047 5.32171763 1.027920e-07
## 3|4 23.4396649994 0.8401688913 27.89875374 2.762585e-171# CIs assuming normality
confint.default(m1)## 2.5 % 97.5 %
## CYTO1 -0.0199697503 2.181642e-01
## CYTO2 -0.0114694136 3.087571e-03
## CYTO3 -0.0040648686 8.732825e-04
## CYTO4 -0.0388016355 1.739791e-03
## CYTO5 -0.0026924641 4.731120e-03
## CYTO6 -0.0010881650 3.294914e-03
## CYTO7 -0.0011874853 3.637983e-03
## CYTO8 -0.0267823980 2.819222e-02
## CYTO9 -0.0186365631 3.495238e-02
## CYTO10 -0.0931519370 4.426732e-02
## CYTO11 -0.0056692013 3.100861e-04
## CYTO12 -0.0703394419 1.745536e-02
## CYTO13 -0.0020599411 9.218450e-04
## CYTO14 -0.0116348580 2.301401e-03
## CYTO15 -0.0766019302 1.220235e-01
## CYTO16 -0.0260478806 5.662181e-01
## CYTO18 -0.0008149684 3.649313e-05
## CYTO19 -0.0293863945 1.104912e-01
## CYTO20 -0.0022347453 6.126248e-04
## CYTO21 -0.0309407550 7.770816e-02
## CYTO22 -0.0053060769 4.445221e-03
## CYTO23 0.0139731765 1.016808e-01#Using step function
step_m1 <- step(m1, trace = F)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredWe also attempt a binary response model with the positive/negative inflammation scores, which seems to be a better fitted model. The confidence intervals for cytokine 4, 9, 10, and 23 do not include 0, which suggests that it may be statistically significant.
xtabs(~ Inflammation, data = dataclean)## Inflammation
## 0 1 2 3
## 27 177 78 18xtabs(~ Necrosis, data = dataclean)## Necrosis
## 0 1 2 3 4
## 75 111 72 36 6# split outcomes (Inflammation or Necrosis Score) with 0-1 vs 2+ as -/+
# try two types of models: ordered factor outcome and binary outcome (try both)# releveled inflammation as factor
dataPO1$Inflammation <- factor(dataPO1$Inflammation)
# binary response model
m3 <- glm(Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO1, family = "binomial")## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(m3)##
## Call:
## glm(formula = Inflammation ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 +
## CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 +
## CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 +
## CYTO20 + CYTO21 + CYTO22 + CYTO23, family = "binomial", data = dataPO1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3855 -0.7277 -0.4042 0.5838 2.5194
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.9594398 0.7981657 -2.455 0.01409 *
## CYTO1 0.1654207 0.0899441 1.839 0.06589 .
## CYTO2 -0.0012708 0.0051856 -0.245 0.80641
## CYTO3 -0.0015651 0.0014822 -1.056 0.29101
## CYTO4 -0.0330060 0.0158206 -2.086 0.03695 *
## CYTO5 0.0069897 0.0120757 0.579 0.56271
## CYTO6 0.0010706 0.0013156 0.814 0.41575
## CYTO7 0.0010141 0.0012318 0.823 0.41033
## CYTO8 0.0313095 0.0292942 1.069 0.28516
## CYTO9 0.0409934 0.0207831 1.972 0.04856 *
## CYTO10 -0.1728190 0.0837291 -2.064 0.03902 *
## CYTO11 -0.0043394 0.0026628 -1.630 0.10317
## CYTO12 -0.0474168 0.0549345 -0.863 0.38805
## CYTO13 -0.0008085 0.0007021 -1.152 0.24951
## CYTO14 -0.0067992 0.0054303 -1.252 0.21054
## CYTO15 -0.1511398 0.1841940 -0.821 0.41190
## CYTO16 0.5637371 0.3519016 1.602 0.10916
## CYTO18 -0.0012869 0.0010949 -1.175 0.23985
## CYTO19 0.0608659 0.0488801 1.245 0.21306
## CYTO20 -0.0020995 0.0025164 -0.834 0.40410
## CYTO21 0.0050385 0.0412417 0.122 0.90276
## CYTO22 0.0039077 0.0033793 1.156 0.24754
## CYTO23 0.0879739 0.0311802 2.821 0.00478 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 125.374 on 99 degrees of freedom
## Residual deviance: 87.069 on 77 degrees of freedom
## AIC: 133.07
##
## Number of Fisher Scoring iterations: 8# CIs assuming normality
confint.default(m3)## 2.5 % 97.5 %
## (Intercept) -3.5238158763 -0.3950636782
## CYTO1 -0.0108665844 0.3417079653
## CYTO2 -0.0114343560 0.0088927965
## CYTO3 -0.0044702025 0.0013400241
## CYTO4 -0.0640138551 -0.0019981919
## CYTO5 -0.0166781949 0.0306575063
## CYTO6 -0.0015078163 0.0036490653
## CYTO7 -0.0014000859 0.0034283666
## CYTO8 -0.0261060666 0.0887251165
## CYTO9 0.0002592515 0.0817276223
## CYTO10 -0.3369249713 -0.0087129802
## CYTO11 -0.0095583507 0.0008795283
## CYTO12 -0.1550864816 0.0602529488
## CYTO13 -0.0021845399 0.0005675757
## CYTO14 -0.0174424859 0.0038440589
## CYTO15 -0.5121534497 0.2098738316
## CYTO16 -0.1259773237 1.2534514546
## CYTO18 -0.0034329049 0.0008590563
## CYTO19 -0.0349374532 0.1566691788
## CYTO20 -0.0070314967 0.0028325405
## CYTO21 -0.0757937711 0.0858707302
## CYTO22 -0.0027156585 0.0105310998
## CYTO23 0.0268619062 0.1490859247We also try a stepwise regression model, which gives the best fittedd model. Cytokines 9, 16, 18, and 23 have p-values less than 0.05, which suggests that they may specifically correlate with immune cell infiltration after transplant.
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftlibrary(leaps)
library(MASS)
# Stepwise regression model
step.model <- stepAIC(m3, direction = "both", trace = FALSE)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred# step.model <- full.model %>% stepAIC(trace = FALSE)
summary(step.model)##
## Call:
## glm(formula = Inflammation ~ CYTO4 + CYTO6 + CYTO9 + CYTO10 +
## CYTO11 + CYTO16 + CYTO18 + CYTO23, family = "binomial", data = dataPO1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5199 -0.8133 -0.5408 0.8247 2.1625
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5454002 0.6126235 -2.523 0.0116 *
## CYTO4 -0.0167166 0.0113717 -1.470 0.1416
## CYTO6 0.0014866 0.0008917 1.667 0.0955 .
## CYTO9 0.0406465 0.0157946 2.573 0.0101 *
## CYTO10 -0.0771683 0.0407839 -1.892 0.0585 .
## CYTO11 -0.0033005 0.0017340 -1.903 0.0570 .
## CYTO16 0.3184029 0.1309310 2.432 0.0150 *
## CYTO18 -0.0013273 0.0005463 -2.429 0.0151 *
## CYTO23 0.0551877 0.0226520 2.436 0.0148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 125.37 on 99 degrees of freedom
## Residual deviance: 101.14 on 91 degrees of freedom
## AIC: 119.14
##
## Number of Fisher Scoring iterations: 7- Dead cells in the liver after transplant? (Necrosis score)
We attempt an ordinal logistic regression model on the numerical necrosis scores. All of the confidence intervals include 0, so none of the cytokines appear to be statistically significant.
library(MASS)
# releveled necrosis as factor
dataPO$Necrosis <- factor(dataPO$Necrosis, levels = c("0", "1", "2", "3", "4"))
dataPO$Necrosis## [1] 0 0 0 0 1 4 2 1 3 1 2 1 1 0 1 4 3 2 3 1 3 2 2 0 1 0 1 2 2 2 2 2 1 3 2 1 1
## [38] 1 1 1 2 1 2 1 1 1 2 3 2 1 0 3 2 1 2 1 3 2 2 0 0 1 1 3 3 2 0 1 1 1 0 0 3 0
## [75] 3 0 0 2 2 0 0 1 1 2 0 0 0 1 1 1 1 1 0 1 0 1 0 0 2 1
## Levels: 0 1 2 3 4# ordinal logistic regression model
m2 <- polr(Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess=TRUE)
summary(m2)## Warning in sqrt(diag(vc)): NaNs produced## Call:
## polr(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 +
## CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 +
## CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO18 + CYTO19 + CYTO20 +
## CYTO21 + CYTO22 + CYTO23, data = dataPO, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## CYTO1 3.192e-02 0.0565505 0.56446
## CYTO2 -2.545e-03 NaN NaN
## CYTO3 1.250e-03 0.0007780 1.60670
## CYTO4 -1.025e-02 0.0098028 -1.04552
## CYTO5 -1.687e-03 0.0016441 -1.02605
## CYTO6 -1.043e-03 0.0009796 -1.06516
## CYTO7 1.455e-03 0.0012184 1.19419
## CYTO8 4.492e-03 NaN NaN
## CYTO9 -2.141e-03 NaN NaN
## CYTO10 -1.346e-02 NaN NaN
## CYTO11 1.323e-03 0.0014214 0.93054
## CYTO12 6.778e-03 0.0146035 0.46417
## CYTO13 -6.015e-04 0.0011060 -0.54383
## CYTO14 -2.169e-03 0.0035172 -0.61663
## CYTO15 -3.605e-02 NaN NaN
## CYTO16 7.691e-02 NaN NaN
## CYTO18 -4.726e-04 NaN NaN
## CYTO19 2.952e-02 NaN NaN
## CYTO20 -1.874e-04 NaN NaN
## CYTO21 3.212e-02 NaN NaN
## CYTO22 4.788e-05 0.0017506 0.02735
## CYTO23 2.922e-02 0.0192031 1.52143
##
## Intercepts:
## Value Std. Error t value
## 0|1 -0.5193 0.5504 -0.9434
## 1|2 1.2675 0.5629 2.2519
## 2|3 2.7644 0.6426 4.3018
## 3|4 4.9783 0.9831 5.0637
##
## Residual Deviance: 259.1091
## AIC: 311.1091## store table
(ctable2 <- coef(summary(m2)))## Warning in sqrt(diag(vc)): NaNs produced## Value Std. Error t value
## CYTO1 3.192063e-02 0.0565505255 0.56446210
## CYTO2 -2.544541e-03 NaN NaN
## CYTO3 1.249995e-03 0.0007779874 1.60670338
## CYTO4 -1.024897e-02 0.0098027847 -1.04551616
## CYTO5 -1.686893e-03 0.0016440617 -1.02605191
## CYTO6 -1.043408e-03 0.0009795748 -1.06516424
## CYTO7 1.454960e-03 0.0012183636 1.19419169
## CYTO8 4.492127e-03 NaN NaN
## CYTO9 -2.140751e-03 NaN NaN
## CYTO10 -1.345868e-02 NaN NaN
## CYTO11 1.322672e-03 0.0014214048 0.93053825
## CYTO12 6.778477e-03 0.0146034506 0.46416955
## CYTO13 -6.014651e-04 0.0011059715 -0.54383420
## CYTO14 -2.168817e-03 0.0035172148 -0.61662903
## CYTO15 -3.604522e-02 NaN NaN
## CYTO16 7.691065e-02 NaN NaN
## CYTO18 -4.725605e-04 NaN NaN
## CYTO19 2.952395e-02 NaN NaN
## CYTO20 -1.874131e-04 NaN NaN
## CYTO21 3.211881e-02 NaN NaN
## CYTO22 4.788159e-05 0.0017505990 0.02735155
## CYTO23 2.921612e-02 0.0192031158 1.52142586
## 0|1 -5.192508e-01 0.5504176604 -0.94337602
## 1|2 1.267495e+00 0.5628545961 2.25190446
## 2|3 2.764412e+00 0.6426239226 4.30175733
## 3|4 4.978260e+00 0.9831229354 5.06372054## calculate and store p values
p2 <- pnorm(abs(ctable2[, "t value"]), lower.tail = FALSE) * 2
## combined table
(ctable2 <- cbind(ctable2, "p value" = p2))## Value Std. Error t value p value
## CYTO1 3.192063e-02 0.0565505255 0.56446210 5.724397e-01
## CYTO2 -2.544541e-03 NaN NaN NaN
## CYTO3 1.249995e-03 0.0007779874 1.60670338 1.081195e-01
## CYTO4 -1.024897e-02 0.0098027847 -1.04551616 2.957845e-01
## CYTO5 -1.686893e-03 0.0016440617 -1.02605191 3.048671e-01
## CYTO6 -1.043408e-03 0.0009795748 -1.06516424 2.868016e-01
## CYTO7 1.454960e-03 0.0012183636 1.19419169 2.324030e-01
## CYTO8 4.492127e-03 NaN NaN NaN
## CYTO9 -2.140751e-03 NaN NaN NaN
## CYTO10 -1.345868e-02 NaN NaN NaN
## CYTO11 1.322672e-03 0.0014214048 0.93053825 3.520925e-01
## CYTO12 6.778477e-03 0.0146034506 0.46416955 6.425263e-01
## CYTO13 -6.014651e-04 0.0011059715 -0.54383420 5.865556e-01
## CYTO14 -2.168817e-03 0.0035172148 -0.61662903 5.374794e-01
## CYTO15 -3.604522e-02 NaN NaN NaN
## CYTO16 7.691065e-02 NaN NaN NaN
## CYTO18 -4.725605e-04 NaN NaN NaN
## CYTO19 2.952395e-02 NaN NaN NaN
## CYTO20 -1.874131e-04 NaN NaN NaN
## CYTO21 3.211881e-02 NaN NaN NaN
## CYTO22 4.788159e-05 0.0017505990 0.02735155 9.781793e-01
## CYTO23 2.921612e-02 0.0192031158 1.52142586 1.281530e-01
## 0|1 -5.192508e-01 0.5504176604 -0.94337602 3.454886e-01
## 1|2 1.267495e+00 0.5628545961 2.25190446 2.432831e-02
## 2|3 2.764412e+00 0.6426239226 4.30175733 1.694488e-05
## 3|4 4.978260e+00 0.9831229354 5.06372054 4.111522e-07# CIs assuming normality
confint.default(m2)## Warning in sqrt(diag(vcov(object))): NaNs produced## 2.5 % 97.5 %
## CYTO1 -0.0789163649 0.1427576217
## CYTO2 NaN NaN
## CYTO3 -0.0002748323 0.0027748222
## CYTO4 -0.0294620747 0.0089641351
## CYTO5 -0.0049091943 0.0015354091
## CYTO6 -0.0029633395 0.0008765233
## CYTO7 -0.0009329891 0.0038429084
## CYTO8 NaN NaN
## CYTO9 NaN NaN
## CYTO10 NaN NaN
## CYTO11 -0.0014632307 0.0041085737
## CYTO12 -0.0218437601 0.0354007144
## CYTO13 -0.0027691295 0.0015661992
## CYTO14 -0.0090624312 0.0047247976
## CYTO15 NaN NaN
## CYTO16 NaN NaN
## CYTO18 NaN NaN
## CYTO19 NaN NaN
## CYTO20 NaN NaN
## CYTO21 NaN NaN
## CYTO22 -0.0033832294 0.0034789926
## CYTO23 -0.0084212983 0.0668535326We also attempt a binary response model on the positive/negative necrosis scores, which gives a better fitted model. Again, all of the confidence intervals include 0, so none of the cytokines appear to be statistically significant.
# releveled necrosis as factor
dataPO1$Necrosis <- factor(dataPO1$Necrosis)
#binary response model
m4 <- glm(Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 + CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 + CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO17 + CYTO18 + CYTO19 + CYTO20 + CYTO21 + CYTO22 + CYTO23, data = dataPO1, family = "binomial")
summary(m4)##
## Call:
## glm(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO3 + CYTO4 + CYTO5 +
## CYTO6 + CYTO7 + CYTO8 + CYTO9 + CYTO10 + CYTO11 + CYTO12 +
## CYTO13 + CYTO14 + CYTO15 + CYTO16 + CYTO17 + CYTO18 + CYTO19 +
## CYTO20 + CYTO21 + CYTO22 + CYTO23, family = "binomial", data = dataPO1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4905 -0.9435 -0.5007 1.0971 2.1463
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0304548 1.2074987 -0.025 0.9799
## CYTO1 0.1051974 0.0757430 1.389 0.1649
## CYTO2 -0.0072987 0.0041645 -1.753 0.0797 .
## CYTO3 0.0009495 0.0013069 0.727 0.4675
## CYTO4 0.0011742 0.0107202 0.110 0.9128
## CYTO5 -0.0124740 0.0208794 -0.597 0.5502
## CYTO6 -0.0014713 0.0013816 -1.065 0.2869
## CYTO7 0.0017876 0.0011727 1.524 0.1274
## CYTO8 -0.0144491 0.0245025 -0.590 0.5554
## CYTO9 0.0097545 0.0163221 0.598 0.5501
## CYTO10 0.0150727 0.0483409 0.312 0.7552
## CYTO11 -0.0001271 0.0018280 -0.070 0.9446
## CYTO12 -0.0138525 0.0240642 -0.576 0.5649
## CYTO13 -0.0008963 0.0006146 -1.458 0.1448
## CYTO14 -0.0085112 0.0047045 -1.809 0.0704 .
## CYTO15 0.0320340 0.0839708 0.381 0.7028
## CYTO16 0.1078550 0.1292101 0.835 0.4039
## CYTO17 -1.0104245 0.7120463 -1.419 0.1559
## CYTO18 -0.0001255 0.0002660 -0.472 0.6372
## CYTO19 0.0400669 0.0378243 1.059 0.2895
## CYTO20 -0.0011347 0.0009620 -1.180 0.2382
## CYTO21 0.0379885 0.0344121 1.104 0.2696
## CYTO22 -0.0007377 0.0027145 -0.272 0.7858
## CYTO23 0.0322030 0.0252215 1.277 0.2017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 132.81 on 99 degrees of freedom
## Residual deviance: 110.52 on 76 degrees of freedom
## AIC: 158.52
##
## Number of Fisher Scoring iterations: 7# CIs assuming normality
confint.default(m4)## 2.5 % 97.5 %
## (Intercept) -2.3971086348 2.3361991026
## CYTO1 -0.0432560885 0.2536508289
## CYTO2 -0.0154610355 0.0008635451
## CYTO3 -0.0016118691 0.0035109007
## CYTO4 -0.0198369429 0.0221854280
## CYTO5 -0.0533967639 0.0284488620
## CYTO6 -0.0041791609 0.0012364694
## CYTO7 -0.0005108973 0.0040861837
## CYTO8 -0.0624731722 0.0335750531
## CYTO9 -0.0222361432 0.0417451927
## CYTO10 -0.0796737313 0.1098191153
## CYTO11 -0.0037099215 0.0034557073
## CYTO12 -0.0610173986 0.0333124625
## CYTO13 -0.0021010175 0.0003083590
## CYTO14 -0.0177319133 0.0007095249
## CYTO15 -0.1325457134 0.1966137291
## CYTO16 -0.1453920737 0.3611021632
## CYTO17 -2.4060095573 0.3851604671
## CYTO18 -0.0006468930 0.0003959766
## CYTO19 -0.0340673446 0.1142011674
## CYTO20 -0.0030202375 0.0007507667
## CYTO21 -0.0294579473 0.1054350298
## CYTO22 -0.0060580481 0.0045825583
## CYTO23 -0.0172302817 0.0816363534We also try a stepwise regression model, which gives the best fitted model. Cytokine 1 has a p-value less than 0.05, which suggests that this cytokine specifically correlates with the dead cells in the liver after transplant.
library(caret)
library(leaps)
library(MASS)
# Stepwise regression model
step.model <- stepAIC(m4, direction = "both",
trace = FALSE)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(step.model)##
## Call:
## glm(formula = Necrosis ~ CYTO1 + CYTO2 + CYTO13 + CYTO17 + CYTO18 +
## CYTO23, family = "binomial", data = dataPO1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5453 -0.9695 -0.6564 1.2369 1.8411
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1873922 0.9237729 -0.203 0.8392
## CYTO1 0.0902847 0.0455846 1.981 0.0476 *
## CYTO2 -0.0041882 0.0024074 -1.740 0.0819 .
## CYTO13 -0.0005261 0.0004327 -1.216 0.2240
## CYTO17 -0.8517493 0.5849126 -1.456 0.1453
## CYTO18 -0.0003877 0.0003022 -1.283 0.1996
## CYTO23 0.0350102 0.0193580 1.809 0.0705 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 132.81 on 99 degrees of freedom
## Residual deviance: 119.89 on 93 degrees of freedom
## AIC: 133.89
##
## Number of Fisher Scoring iterations: 6